%% Actual estimation. Load data
clear
% workdir = '/Applications/Dropbox/RnD_KK/tradecost_replication';
workdir = '/Dropbox/RnD_KK/tradecost_replication';

% addpath = strcat(workdir,'/matlab_model2');
path(workdir,path);
addpath(genpath(strcat(workdir,'/ReStat_code')))

prm.s = 6;
thresh = '20';               % # entrants in a destinanation must be >= thresh
data.year = '2004';

M = csvread(strcat(workdir,'/data_matlab/data_qty_',thresh,'thresh_',data.year,'.raw'));           % (IxP)xK
M2 = csvread(strcat(workdir,'/data_matlab/data_price_',thresh,'thresh_',data.year,'.raw'));           % (IxP)xK

% Load number of firms/obs per product
% numfirms is a vector with [product_id number of firms sigma]
data.numfirms = csvread(strcat(workdir,'/data_matlab/data_',data.year,'_count_',thresh,'thresh.raw'));  % Px1

data.P = size(data.numfirms,1);
data.K = size(M,2)-2;           % # countries
K = data.K; params.K = data.K;

% Data is originally in "long" format, so products appear sequentially.
% Rearrange to "wide" format.
% PS: We need to make sure that all products in the CSV file are actually
% used in estimation. I.e. that no products have less than $thresh firms in
% all destinations.

thresh2 = str2num(thresh);
start = 1;
for p=1:data.P
  stop = start+data.numfirms(p,2)-1;
  qty = M(start:stop,3:end);    % I_k x K Sales (quantity)
  price = M2(start:stop,3:end);
  entrants  = qty>0;            % (I_k x P) x K Export participation
  
  for n=1:data.K
    tmp = qty(:,n);
    tmp_p = price(:,n);
    if length(tmp(entrants(:,n)))>=thresh2           % Check if we want to use this product-dest pair
      data.usedata(n,p)=1;
      data.sales_actual{n,p} = tmp(entrants(:,n));
      data.lnX{n,p} = log(data.sales_actual{n,p});
      data.lnX_dev{n,p} = data.lnX{n,p} - mean(data.lnX{n,p});

      data.price{n,p} = tmp_p(entrants(:,n));
      
      data.priceinv{n,p} = 1./data.price{n,p};
      data.priceinv_dev{n,p} = data.priceinv{n,p} - mean(data.priceinv{n,p});      

      data.priceinv2{n,p} = data.price{n,p}.^(-2);
      data.priceinv2_dev{n,p} = data.priceinv2{n,p} - mean(data.priceinv2{n,p});      

      data.priceinv3{n,p} = data.price{n,p}.^(-3);
      data.priceinv3_dev{n,p} = data.priceinv3{n,p} - mean(data.priceinv3{n,p});      

      data.priceinv4{n,p} = data.price{n,p}.^(-4);
      data.priceinv4_dev{n,p} = data.priceinv4{n,p} - mean(data.priceinv4{n,p});      

      data.lnprice{n,p} = log(data.price{n,p});
      data.lnprice2{n,p} = data.lnprice{n,p}.^2;
      data.lnprice3{n,p} = data.lnprice{n,p}.^3;
      data.lnprice4{n,p} = data.lnprice{n,p}.^4;

      data.lnprice_dev{n,p} = data.lnprice{n,p} - mean(data.lnprice{n,p}); 
      data.lnprice2_dev{n,p} = data.lnprice2{n,p} - mean(data.lnprice2{n,p}); 
      data.lnprice3_dev{n,p} = data.lnprice3{n,p} - mean(data.lnprice3{n,p}); 
      data.lnprice4_dev{n,p} = data.lnprice4{n,p} - mean(data.lnprice4{n,p}); 
      
      data.priceinvlnp{n,p} = data.priceinv{n,p}.*data.lnprice{n,p};
      data.priceinvlnp_dev{n,p} = data.priceinvlnp{n,p} - mean(data.priceinvlnp{n,p});      
      
      data.lnValue{n,p} = data.lnX{n,p} + data.lnprice{n,p};
      data.price_mean(n,p) = mean(data.price{n,p});
      data.price_median(n,p) = median(data.price{n,p});
      
      % Calculate phi
      tmp = cov(data.lnX{n,p},data.lnprice{n,p});
      data.var_lnX(n,p) = tmp(1,1);
      data.var_lnprice(n,p) = tmp(2,2);
      data.cov_px(n,p) = tmp(1,2);
      data.corr_px(n,p) = corr(data.lnX{n,p},data.lnprice{n,p});

      phi(n,p) = tmp(1,2)./tmp(2,2) + prm.s;
      phisig(n,p) = tmp(1,2)./tmp(2,2);
      
      [data.ECDF{n,p},data.esales{n,p}] = ecdf(data.sales_actual{n,p});    
      data.entrants(n,p) = length(data.sales_actual{n,p});
      
      % Construct pct ratios
      s = data.lnX{n,p};
      pcts = prctile(s,[5:5:95]);
      if (pcts(1)==pcts(2)) pcts(2)=pcts(2)+.1; end;
      if (pcts(2)==pcts(3)) pcts(3)=pcts(3)+.1; end;
         
      data.emom(n,p,1) = (pcts(end)-pcts(end-1)) / (pcts(2)-pcts(1));
      data.emom(n,p,2) = (pcts(end-1)-pcts(end-2)) / (pcts(3)-pcts(2));
      data.emom(n,p,3) = var(s);

    else data.usedata(n,p)=0; end;
  end  
  start = stop+1;
end

% phi should only be product specific
data.phi = mean(phi);        % 1 x P
data.phisig = phisig;        % N x P   

% Number of destinations per product
data.K2 = sum(data.usedata)' 
data.avgentrants = mean(data.entrants(logical(data.usedata)));

clear M N;
disp('Year'); disp(data.year);
disp('Number of products'); disp(data.P);
disp('Number of countries'); disp(data.K);

disp('Number of products usable in estimation per destination');
disp(sum(data.usedata'));

% Number of usable product-destinations
data.num_proddest = sum(sum(data.usedata));
disp('Number of product-destinations'); disp(data.num_proddest);

if any(sum(data.usedata)==0) error('Some products are have no useable destination'); end;
disp('Data loaded');

%% NLS Joint estimation

res = []; lb = []; ub = [];
K = data.K;
P = data.P;

%
% Constraints and initial values
%

%tmax = log(data.price_median(data.usedata==1)'*10);

prm.idx_t_n = 1:data.K;
prm.idx_t_p = data.K+1:data.K+data.P-1;
prm.idx_phisig1_p = data.K+data.P:data.K+2*data.P-1;
prm.idx_a_p = data.K+2*data.P:data.K+3*data.P-2;
prm.idx_a_n =data.K+3*data.P-1:2*data.K+3*data.P-2;

b0(prm.idx_t_n) = 2;
b0(prm.idx_t_p) = 0;
b0(prm.idx_phisig1_p) = -2;
b0(prm.idx_a_p) = 1;
b0(prm.idx_a_n) = 1;

lb(prm.idx_t_n) = -20;
lb(prm.idx_t_p) = -20;           
lb(prm.idx_phisig1_p) = -20;     
lb(prm.idx_a_p) = -20;     
lb(prm.idx_a_n) = -20;

ub(prm.idx_t_n) = 20;           
ub(prm.idx_t_p) = 20;
ub(prm.idx_phisig1_p) = 0;       
ub(prm.idx_a_p) = 20;     
ub(prm.idx_a_n) = 20;


% Linear inequality: lnt_k+lnt_n  < ln(median price) & -lnt_k-lnt_n < -ln(0.01*median price)
tmp = data.price_median(:,1:data.P-1);
num_kn = sum(tmp(:)>0);
AA = zeros(num_kn,length(b0));
bb = zeros(num_kn,1);
ineq = 1;
for kk = 1:data.K
   for pp = 1:data.P
     if data.usedata(kk,pp)==1
       AA(ineq,prm.idx_t_n(kk)) = 1;
       if pp<data.P 
           AA(ineq,prm.idx_t_p(pp)) = 1;
       end     
       bb(ineq,1) = log(1*data.price_median(kk,pp));      
       ineq = ineq+1;
     end
   end
end

ineq = ineq+1;
for kk = 1:data.K
  for pp = 1:data.P
    if data.usedata(kk,pp)==1
      AA(ineq,prm.idx_t_n(kk)) = -1;
      if pp<data.P
        AA(ineq,prm.idx_t_p(pp)) = -1;
      end
      bb(ineq,1) = -log(.01*data.price_median(kk,pp));
      ineq = ineq+1;
    end
  end
end



f = @(b) imo_2011_NLS_all3FE(b,data,prm,0);

options = optimset('Display','iter','TolFun',1e-8,'TolX',1e-8,'MaxFunEvals',5e7,'MaxIter',10000);
%options = optimset(options,'Algorithm','interior-point');
options = optimset(options,'Algorithm','sqp');

tic, [res.b,res.fval,res.exitflag,res.output,res.lambda,res.grad,res.hessian]  = fmincon(f,b0,AA,bb,[],[],lb,ub,[],options); toc

res.lnt_n = res.b(prm.idx_t_n);
res.lnt_p = [res.b(prm.idx_t_p) 0];
res.phisig1_p = [res.b(prm.idx_phisig1_p)];
res.a_p = [res.b(prm.idx_a_p) 1];
res.a_n = res.b(prm.idx_a_n);

res.lntmat = repmat(res.lnt_n',1,data.P) + repmat(res.lnt_p,data.K,1);
res.a = repmat(res.a_n',1,data.P).*repmat(res.a_p,data.K,1);
sderr = sqrt(diag(inv(res.hessian)))';

res.lnt_n_sd = sderr(prm.idx_t_n);
res.lnt_p_sd = [sderr(prm.idx_t_p) 0];
res.phisig1_p_sd = [sderr(prm.idx_phisig1_p)];
res.a_n_sd = sderr(prm.idx_a_n);
res.a_p_sd = [sderr(prm.idx_a_p) 0];

disp('lnt_n'); disp([res.lnt_n' res.lnt_n_sd']);
disp('lnt_p'); disp([res.lnt_p' res.lnt_p_sd']);
disp('phisig1_p'); disp([res.phisig1_p' res.phisig1_p_sd']);

disp('Estimated trade cost relative to median price');
disp(exp(res.lntmat)./data.price_median);

%save(strcat(workdir,'/matlab_model2/res_2004_20thresh_Restat7'),'res'); 
